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An overview of some recent developments in inhomogeneous models is presented. 
As the volume and precision of cosmological data improves, it will become more and more essen- 
tial to understand the non-linear behaviour of the Einstein field equations. This requires the study 
of exact inhomogeneous solutions, including their density distributions, their evolution, their ge- 
ometry, and their causal structure. Observations are strongly affected by the detailed geometry 
and evolution of a model, and therefore interpretation of observations depends on understanding 
them. 

It is generally assumed the universe is homogeneous if averaged over large enough scales, but to 
actually prove this is so, will require the assumption to be relaxed, and a rigorous inhomogeneous 
approach to be applied. 

Though the Lemaitre-Tolman metric has long been used for models of spherical inhomogeneities, 
there have been a number of new results, including a variety of methods for creating models with 
specific properties, and their apphcation to cosmic structures on several different scales. 
Interest in the Szekeres metrics is on the increase, and the quasi-spherical metric was recently used 
to model specific cosmic structures for the first time. The quasi-planar and quasi-hyperspherical 
metrics have been hardly studied until recent work invesigated their physical and geometric prop- 
erties. There is enormous scope for work with these metrics. 



5th International School on Field Theory and Gravitation, 
April 20-24, 2009 
Cuiaba city, Brazil 



© Copyright owned by the author(s) under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike Licence. 



http://pos.sissa.it/ 



Modelling Inhomogeneity in the Universe 



1. Introduction 



Why study inhomogeneous models? The real universe is very lumpy. To properly understand 
what we see, we should apply all possible methods: perturbation theory, A^-body newtonian sim- 
ulations, and exact inhomogeneous metrics — each has its domain of validity. Inhomogeneous 
metrics have the advantage that they are fully non-linear and relativistic solutions of the Einstein 
field equations (EFEs). 

The assumption of homogeneity has become so well established, that it has become all- 
pervasive. But now, with so much data coming in, it's time to test homogeneity. Cosmological 
data reduction relies heavily on the Robertson-Walker (RW) metric — we need to beware of a cir- 
cular argument. It will be a significant challenge to check which of our well-known results actually 
depend on the assumption of homogeneity, and to re-derive them all without that assumption. 

Here I will present a selection of results in inhomogeneous cosmology, especially work done 
with Lu, McClure, Krasiriski, Bolejko, Celerier, Alfedeel, Mustapha, Ellis and others, but I won't 
try to be comprehensive. I'll attempt to provide the basics, and thereby promote the use of inho- 
mogeneous metrics for the study of cosmological problems. 

Inhomogeneous metrics will become more important as the amount and accuracy of cosmo- 
logical data increases, and more precise analysis is needed, so there are plenty of opportunities for 
good research. 



2. The Lemaitre-Tolman Metric 



The Lemaitre-Tolman (LT) metric was the first inhomogeneous non-vacuum metric to be dis- 
covered, and has probably been the most popular choice for modelling cosmic inhomogeneity ever 
since, certainly in recent decades. It is a spherically-symmetric, inhomogeneous dust model, dis- 
covered by Lemaitre, rediscovered by Tolman, and studied by Bondi [84, 115, ^]. The metric 
is 



-dr + 



1 + / 



(2.1) 



where dfJ^ = dO"^ + sin^ 6 dc/)'^, R{t, r) is the areal radius, and R' = dR/dr. The free function /(r) 
determines the local geometry; it gives the "embedding angle" of constant t,9 surfaces in 3-d flat 



space [pl|]. Also the Ricci scalar of the spatial 3-surfaces, 

^_ -2{Rf + fR') 

Rm' 

is only zero (for all r) if both / and /' are zero. The matter is a pressure-free perfect fluid, 
that is comoving 



(2.2) 



(2.3) 



(2.4) 
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From the EFEs we get 



where R = dR/dt, and 



r:' 



2M , Ai?2 



2M' 



Kp - 



(2.5) 



(2.6) 



R'^R' ' 

where M(r) is a second free function that gives the gravitational mass within a comoving shell of 
radius r. Here /(r) also plays the role of twice the local energy per unit mass of the dust particles, 
so it's often written /(r) = 2E{r). It follows from (2.5) that 



■■ M AR 



When A = 0, the solutions of (^), in terms of parameter r], are 
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A = 0, / >0: 
A = 0, / = 0: 
A = 0, /<0: 
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R 



R 
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(coshry — 1) 
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(2.7) 
(2.8) 

(2.9) 
(2.10) 
(2.11) 



for hyperbolic, parabolic, and elhptic evolution respectively. (Near the origin, where / — > 0, the 
type of evolution is determined by the sign of Rf /M or //M^/^.) When A / there is a very com- 
plicated solution in terms of elliptic integrals. These solutions contain a third free function a(r), 
which is the local time of the big bang, the time when i? = on each worldline. In other words, 
the constant r worldlines all emerge from the bang at different times, usually the outer spheres first 
and the origin last, as illustrated in the sketches below. 




The above evolutions equation and solutions may also be written 

■ 2 Aa2 

/3 = - + x + — - , 
a 6 



(2.12) 
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A = 0, / >0: 

A = 0, / = 0: 

A = 0, / < 0, < 7? < vr : 

A = 0, / < 0, TT < 7? < 27r : 

where 



t = a + I v'(l + xa)2-l - arcosh(l + xa) | , (2.13) 
^/2a3 

} , (2.15) 



a + - 



t = a + x 



-3/2 



arccos 



l + xa)- Vl-(1 + 



t = a + x ^/^|27r-arccos(l + xa) + (l + xa)2 | . 

(2.16) 

R \f\ ^ R 



a 



(2.17) 



Naturally, the time reverses of these models, obtained by changing {t — a) to (a — t) and a' to 
—a', are also solutions. It is quite possible to have adjacent elliptic and hyperbolic regions in 
one model — for example, a re-collapsing dust cloud could be surrounded by an ever-expanding 
universe. The two regions would have a parabolic shell at the boundary between them, but extended 
paraboUc regions are also possible. In practice. 



( P.IID , ( I2.13D , ( I2.15D , ( p.l^ ) are not good for 
calculating the evolution of worldlines that are close to parabolic, so a series expansion is used 
instead. Similarly, near the bang or crunch, where the evolution is close to parabolic, one obtains 
better accuracy by using the same series expansion. 



It is sometimes useful to have an expression for R'. When A = it follows from (p^-(|2!TT|), 



that for all / values one can write [|63|] 



\M f 

Alternatively, one can write the parametric expressions 



, , M' 3/' , , 
a + —](t-a) 
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(1-01) + ^(^01 
sin 77(77 — sin 77) 



M 



(1- 



(1 — COS7/)2 



sm 



(1 — cos 77)^ 



sinh77(sinhr/ — 77) 



(cosh 77 — 1)2 
A scale length and time may be defined by 

M 



Rir) 



f{r) 



4>5 



M 



sinh77 
(cosh 77 — 1) 



(2.18) 

(2.19) 
(2.20) 
(2.21) 
(2.22) 

(2.23) 



and for elliptic worldlines the maximum R is 2R, while the lifetime from bang to crunch is 27rT. 

By specifying A and the three free functions — M{r), f{r), and a{r) — an LT model is fully 
determined. Between them they provide a radial co-ordinate freedom and two physical relation- 
ships, e.g. M = M{r), f = f{M) and a = a{M), though it is normal to give all of them in terms 
of r. It is not possible to give any kind of standard form for one of these functions that will cover 
all possibilities. For example, the choice M oc r'^ is common, but does not allow regions of vacuum 
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where M' = 0; a standard choice for /(r) cannot include both models in which / changes sign, 
and those in which it doesn't; and similarly no choice of a(r) can cover cases where a is constant 
in some places and cases where it never is. 



See [ ]74[ ] for a survey of work done on inhomogeneous models up to 1997, QlOll ] for an intro 



duction to some inhomogeneous models, and [20] for a summary of some recent developments. A 



dynamical systems analysis is given in [109] 



2.1 Singularities 

Singularities occur where the density (|2.6|) or the curvature diverge. The Kretschmann scalar 



IS 



^ P rjabcd 48M2 32MM' 12{M') 



(2.24) 



Big Bang At the big bang or the big crunch, we have R = 0, which occurs where t = a or where 
t = a + 2'itT. Here R' diverges unless a' = 0. The bang and crunch surfaces are spacelike [^, |6^ ], 
except possibly at the origin. 

Shiell Crossings Shell crossings are timelike surfaces that occur where an inner spherical 
shell of matter collides with an adjacent outer shell, so that R' = 0. These surfaces are timelike 
[50, and have a different redshift structure from the bang [|63|]. Since the r coordinate is 
comoving, it becomes degenerate at such loci. Physically one might argue that non-zero pressure 
would develop before a shell crossing occurs, but for a "fluid" of many stars or galaxies that doesn't 
apply. Clearly shell crossings represent a breakdown of the LT assumptions and for many purposes 
they are undesireable. Shell crossings can be eliminated from an entire model, in the A = case, 
by applying the conditions found in [|^, to the 3 arbitrary functions. These conditions were 
derived by writing R' in terms of the parameter rj and looking at the early and late time behaviours. 
They are important if you want your model to be everywhere well behaved.^ 



rlier 





■gy not always increasing outwards 
and R' = can occur at non-singular locations as explained below. 



However, both R 

'Shiell Focussing' There are also "shell focussing" singularities, e.g. [I37|, B2 



8 102, 65 



|83| , |49| , 81], 72]. For certain LT models, the first event of the big crunch to form, where the 
central worldline reaches the crunch surface, can emit many light rays, some of which may even 
reach infinity. (So they might be better called "light focussing" singularities.) The nature of the 
singularity is difficult to understand, and seems to depend on the path of approach to the singular 
point. 

12lt|97|,Pl 



Shell crossings have been extensively investigated, e.g. 
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2.2 Regularity Conditions 



Regular signature For the metric (gJJ) to retain a Lorentzian signature, 

/>-l (2.25) 

is required, the equality only occuring where R' = — see below. 
Reguiar Origins An origin of spherical coordinates is a locus Tq where 

R{t,ro)=0 V t , (2.26) 

so that R{t,ro) = 0, R{t,ro) = 0, etc. Obviously one usually wants an origin to be a normal 
timelike worldline. The conditions for a regular origin are obtained by requiring that, in the limit as 
the origin is approached, the density and the curvature should not diverge, and the time evolution 
at the origin should be a smooth continuation of it's immediate neighbourhood. See for example 



8|, |94|]. It is found that, away from the bang or crunch, on a constant t slice, 

Mr^R^ , f^R^. (2.27) 
This may be realised by setting / oc M^/^, e.g. M ^ r^, f ^ r^, at the origin. Variables a, x & /? of 



(2.17) have the advantage that they are non-zero at the origin. If in addition one wants the density 
to be smooth through the origin, i.e. to have zero gradient there, then there are further conditions 
[f4|l> most notably 

a' ^ . (2.28) 

However, there is no singularity if this last one does not hold. Thus, the locus R = includes both 
the spaceltke bang and crunch surfaces, and the timelike origins. 



Reguiar Spatiai Extrema Similarly, R' = includes regular loci as well as singular shell 



crossings. As pointed out in [122], any spherically symmetric model with closed sptial (t = const.) 
sections, such as the A; = +1 Freidmann-Lemaitre-Robertson-Walker (FLRW) model, has an ori- 
gin, a maximum radius, and a second origin — a north pole, an equator and a south pole. At a 
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spatial maximum we obviously have R' = 0, but we expect the density and curvature to be regular 
at such a locus. This is only possible on a comoving shell, i.e. 



V i . 



(2.29) 



The conditions for a regular maximum [64, 26] are that there is no shell crossing and no surface 
layer, i.e. 



M'{rra) = = firm) = a\rra) , firm) 



-1 . 



(2.30) 



Therefore the LT model may have a number of interesting spatial topologies \ fi4[ pl| , [78| ], such as 
a black hole in a cosmological background, or a sequence of maxima and minima — "bellies" and 
"necks". It is also possible to have an elliptic (recoUapsing) model that is open. 

2.3 Special Cases 

Dust Robertson-Walker The LT metric contains the dust RW metric as the special case 

/ oc m2/3 , a' = 0. (2.31 ) 

Putting this in ( |2.5[ ) and (2^) makes R/M^^^ and p independent of r. In standard RW coordinates, 
M = (Kpo'S'o/6)r^, / = —kr^, a = 0, R = rS'(t), so it is evident that Sit) is the scale factor, and 



/2 



-k 



6 



(2.32) 



Consequently one may write the LT arbitrary functions in a form that looks like RW plus perturba- 
tion, but is exact. 



M = Msr^il + Mir)) 
f = f2r\l + fir)) 
a = ao(l + a(r)) , 



(2.33) 
(2.34) 
(2.35) 



where M, / and d may be set to zero at the origin, say. In terms of the RW parameters of the 
'unperturbed' RW model that applies at the origin, we can write 



S — — , H — — , Qr. 

r b 

so that f2 = —k — > 

and of course Sq 



^ 2M3 

S3H2 

M3 

1 



f2 



52^2 



A 

3IP 



mO 



2i/2(_A,f],o)3/2 



ao = , 



HoV—kTlko 



A = mAH^ 



(2.36) 
(2.37) 

(2.38) 



Schwarzschlld The spherical vacuum metric is obtained if M' = 0, and the different choices 
of fir) and a(r) cover it with different families of geodesic coordinates. But to get the full 
Scwarzschild-Kruskal-Szekeres (SKS) topology requires / = — 1 and a' = = /' so that i?' = at 
the "throat" or "neck", and that a decreases and / increases on either side — see [51, 54] for the 
details and some plots. 
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Asymptotic 
■< — region 



f(r)/ 




Asymptotic 
region — >• 




Datt-KantOWSki-Sachs The Datt [ p5[ ] models are inhomogeneous Kantowski-Sachs type mod- 
els, and though often treated as a separate solution with R' = 0, they are in fact limits of LT models 

111- 

Vaidya In the null limit, when / — > cx) we get the Vaidya metric that represents incoherent 



radiation emanating from (or converging on) a spherical body [pq, 53] 



2.4 Constructing Inhomogeneous Models 

The most obvious way to construct an LT model is to choose the three arbitrary functions. 
Choosing /(r), for example, works quite well if one is interested in the geometry and topology 



of the model. See [|64|, qJJ]. But for many situations, it is not always obvious what the density 
distribution and evolution will be, given /, M and a. In [ |1 10| ], for example, the use of the density 
Pi{r), the 3-d Ricci scalar ^i{r) and the areal radius Ri{r), on an initial surface at t = ti was 
advocated, and an appendix suggested how 'lumps' and 'voids' in the density and curvature could 
be prescribed on the initial surface. 

In place of a(r), one may instead specify Ri = R{ti,r) at some initial time, set = at t = t j 
and re- write ( 2.1 1[ ) in the form 



R 



M 

FT) 

M 



(1 — COST]) + Ri COST/+ , 



2M 



1 sin 7] 



[t] — sinrj) + 



Ri 



sinr/+ . 



2M 
-f)Ri 



1(1 — COS?]) 



(2.39) 



etc for the other cases, (2.10) & (2.9) 



If you prefer to think in terms of R, then choose M{Ri), f{Ri) and a(i?i) on an initial surface 
t = ti, and set r = iij. If p = pi{Ri) is given, then, again choosing r = Ri, 



M{Ri)-Mo 



Hcp,{R)R' 



dR. 



(2.40) 



Ro 
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If the expansion rate and radius Ri{M) and Ri{M) are specified, then by (2.5) 



Ui 



AR? 



(2.41) 



or if the Ri{M) and f{M) are specified, then 
«(M), + 

or R{M) : 



AX 



X= ( V9M2A-(i?2-/)3 - V9M2A j 



\ 1/3 



2M 



when A = . 



(2.42) 
(2.43) 



Initial and Final Profiles In ^ ^ a number of very useful methods for A = LT models 
were presented. Since these have been summarised elsewhere 1 20 ] we will only outline the basic 
idea here. Rather than specifying all the data on a single 'initial' 3-surface, one may instead specify 
the density profile pi on one constant time surface t = ti and another density profile p2 at a later 
time t = t2. There is a well-defined algorithm for finding the LT model that evolves from one to 
the other. 

Suppose, on the surfaces t = & t = we specify the density to be p = pi{M) and p = 
P2(M),then2from(0), 



fM g 

RKM)-RI= —TTTsdM, 



1,2 



(2.44) 



and normally we would have Rq = = Mq. We set the coordinate freedom via r = M. Since 
M is constant along each particle worldline, we now know Ri and i?2 for each particle. We 
consider a specific M, and we assume R{ti,M) > and R2> Ri- By the time t2, the worldline is 
either hyperbolic and still expanding (HX), elliptic and still expanding (EX), or elliptic and already 
collapsing (EC). In the HX case, we apply (2.13) and (2.17) at the two times and subtract them: 



+ xa2)'^ — 1 — arcosh(l + xa2) 



+ xai)2 — 1 + arcosh(l + xai) — x ^^'^{t2 — ti) = i^Hxix) = , 



(2.45) 



with similar expressions for the other cases. This is solved numerically by the bisection method, 
and for this purpose, a pair of x values that bracket the solution were found. Having obtained 
/ = xM2/3, a is found by using ( |2.13| ) again: 



ti-x ^/^|\/(l + xai)2-l -arcosh(l + xai)| . 



(2.46) 



Obviously it is important to know which case applies along each worldline. In [ J76| ] it was shown 
that 



t2-tl> 



02 \ 

y) 



vr — arccos 1 



2ai 

a2 



+ 2i 



02 



EC 



(2.47) 



In this case, though, ( 2.44 ) would not be well defined if there were vacuum pi{M) anywhere, since the range of 7? 
over which M is constant could be anything. 
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t2-ti<^(al^^-a/^) HX (2.48) 

otherwise EX (2.49) 



The borderUnes between these cases require careful treatment, see appendix B of []61|]. If R2 < Ri 



then there are 3 more solutions, including the collapsing hyperbolic model — the time reverse of 



(2.13). 



A similar approach may be used if velocity profiles Ri{M) and i?2(M) are given at ti & t2, 
of if a density profile is given at one time and a velocity profile at another [[77|]. There are quite a 
few other useful options, such as setting the late time density of velocity behaviour, specifying a 
simultaneous time of maximum expansion, specifying only growing or decaying modes, etc [jST]]. 

Applications of these methods to model a galaxy cluster, a void, a galaxy with a central black 



hole, the Shapley Concentration and the Great Attractor, etc, can be found in [ jTq , 19,13,1^]. 

3. Observations in Lemaitre-Tolman Models 

The assumption that the universe is homogeneous, and thus well-represented by an FLRW 
model, has led to a very good understanding of its very large scale features and evolution. But once 
the cosmological data are sufficiently accurate and complete over a large enough range of redshifts, 
this assumption should be checked. 

However, the assumption of homogeneity pervades so much theoretical and observational 
work so thoroughly, that there is a real danger of a circular argument. Consequently, any proof 
of homogeneity must ensure it does not rely on results obtained using an assumption of homogene- 
ity. Clearly this will not be a simple task. More precisely, the aim is not only to verify homogeneity 
but also to quantify it: how much fluctuation is there on each scale? 

There are a several reasons why spherical symmetry is a good first step towards relaxing the 
homogeneity assumption: (a) we are at the centre of our past null cone, so it makes sense to consider 
spacetime in terms of spherical co-ordinates about the observer; (b) the universe does seem close 
to isotropic on large scales, but radial homogeneity is not easy to verify because of the finite travel 
time of light and the miniscule duration over which cosmological observations have been made, so 
it is more urgent to determine the radial variation of the metric; (c) there is no deep all-sky redshift 
survey at present, and the zone of avoidance is likely to be a gap in any survey for the foreseeable 
future; (d) it keeps the theory and numerics tractable while the basics are sorted out. Of course, in 
the long run, the assumption of spherical symmetry will be dropped. 

Here we derive the observational relations that would be expected in an LT model with given 
arbitrary functions. This is also known as the 'forward problem'. Below we focus on a central 
observer, though non-central observers have also been considered. 

3.1 Observables and Source Evolution 

The observables we shall use are those for which the dataset is already substantial and will 
in the near future become extensive, the redshift z, the number density in redshift space n{z), the 
apparent luminosity and angular diameter £{z) & 5{z). Connected with each of these is a source 
property, the peculiar velocity C,, the mass per source jjL{z) the absolute luminosity L{z) and the 
true diameter Diz). 
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The redshift z is 



- , 1 (3.1) 

where Ao and Ag are the observed and emitted wavelengths. The diameter and luminosity"^ distances 
are 

dD = j, dL = ^dio = W^^-'^y'dio , (3.2) 

where 6 and £ are the angular diameter and apparent luminosity of a source, D{z) and L{z) are 
the corresponding true diameter and absolute luminosity, m and m are the apparent and absolute 
magnitudes, and diQ is 10 parsecs. The two distances are related by the reciprocity theorem [44, 

Toq,^, 

{l + zfdD = dL. (3.3) 



The Hubble and deceleration constants are obtained from the slope and concavity of the diiz) plot 
at the origin. 



ddi 



dz 
1 d^dL 



1 



Hq dz2 

and a common observational definition of H{z) and q{z), based on the FLRW model, is 

1 1 d f di 



(3.4) 
(3.5) 

(3.6) 



H ^l + Vtk{HodL/{l + z)Y dz\{l + z 

For general non-homogeneous models, there is no obvious general definition of H{z) or q{z), and 
a number have been proposed. In any case, what matters is the relation between the model and 
observations. 

If in a redshift survey of the sky, dN sources are observed to lie between z and z + dz within 
solid angle du = sin 9 dO dcj), then the redshift-space mass density is 

2a fidN 
K dujdz 

where n is the redshift space number density and fi is the mean mass per source. For a treatment 
which considers several different source types and observations at different wavelengths see [|^. 

A significant feature of these definitions is that each observable, 6, i and n, is associated with a 
source property, D, L and fi, which have certainly evolved over cosmological timescales. The latter 
are much harder to determine observationally, and studies of their values and evolution invariably 
assume a homogeneous RW model in which to do the analysis. However, if we eventually want to 
prove that the universe is homogeneous, it is imperative to avoid a circular argument. The only way 
to be certain of the conclusion is to do the analysis without making the homogeneity assumption. 

■^In Jsc] ] a "corrected luminosity distance" was defined to be the same as the diameter distance. Some authors 
have called this latter the "luminosity distance", which perhaps has led to a confusion of terminology and sometimes to 
incorrect definitions. 
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3.2 The Null Cone and the Observational Relations 

Light rays arriving at the central observer O follow ds'^ = = d9'^ 
past null cone (PNC) of the observation event (t = tQ,r = 0) satisfies 



dr 



w ' 



W = Vl + 2E , 



% so from (2.1) the 



(3.9) 



and we write the solution t = i{r) or r = r(t), defining the local time from the bang to O's PNC 
with 



T = t — a . 



(3.10) 



This radial null path is necessarily geodesic. We denote a quantity evaluated on the observer's past 
null cone with a hat on top or as a subscript, for example R{i{r),r) = Roi [R]/^, though this will 
often be omitted where it is obvious from the context. For a given LT model, equation (3^) must 
be solved numerically. 

As is well known for the LT model (e.g. [23, |95|, 0), the redshift of sources on the PNC 
observed at O obeys 



dz 



(1 + ^) 



R/_ 



dr 



(3.11) 



where R' is given by (2.8) and (2.5). 



The diameter distance is, by ([3.2|), the quantity that converts measured angular sizes of objects 
to their physical sizes at the time of emission. It is evident from the metric ( p] ) that this is the 
areal radius R, evaluated on the PNC, 



dD = R = R{i{r),r) . 



(3.12) 



and of course di follows from (3.3). 

To convert the proper density of an LT model to the observed redshift space density, requires 
that we know how the locus of the PNC relates z to comoving radius, i.e. f{z). Then the total mass 
contained in a small volume must be the same: 



— dz duj 

K 



/7?2 



pR'R 
W 



■ drdio 



-2 dz 
npR =2a— and a 
dr 



M' dr 
Wd^ 



(3.13) 



where ( |2.6[ ) has been used. 

In the numerical solution of these equations, (3^), ( 3.11 ) & ( 3.13 ), we need to evaluate R, R' 
and R! at each new point along the PNC. Along the constant r worldline at each step we integrate 



and 



^dR 
~R 



dR' _ RR! 
~dR~^ 



(3.14) 



(3.15) 
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where R, and RR' come from and ( |2.8[ ). Equations ( p.l4| ) and ( p.l5[ ) are solved in one 
numerical integration for each step of integrating (3^), ( 3.11 ) & ( 3.13 ).^ 



The variation of the areal radius down the PNC is 



= 

dr dr 



R' \l 



R 
W 



(3.16) 



and its second derivative is 

d^^ f d di d\dR 

+ ■ 



dr^ \ dr dr dt J dr 
R" 



W 



, R\ R' RW RR' 
1 H -i? H \ 

w w\ w w 



(3.17) 
(3.18) 



w ] w 3 Ry \ w \ 



M'R'R 
W^R 



(3^9) 



where W' and R were eliminated using (2^) and {21). It is important for later to note that R{r) 
may have a maximum value where di2/dr = 0, and at this locus we have 



R = W ^ 



2M , , , 



and consequently, using (3.13), 



d?R aR' dz 



RW dr 



Thus the slope of the diiz) curve is 

ddL d{l + zfR 



dz 



dz 



2{l + z)R + {l + z 



R' 



^2dR/dr 
dz/dr 



and at the origin we have 



{l + z)\2R + —{W-R)j , 



1 _ R' 



R' 



z=0 



The definitions for the radial and tangential Hubble rates 



Hr — , Hf — — 

R' ' R 



The A = special case is much easier, because (2.18) allows one to integrate (3.9) without solving (3.14) at every 



(3.20) 



(3.21) 



(3.22) 



(3.23) 



(3.24) 



step. 
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represent the metric expansion rates in the radial and tangential directions, but one needs to be 
careful which of these, if any, relates to which observation. 

Near the origin — the vertex of the PNC — in addtion to ( 2.26 ) & ( 1.21 ), we have 



W 







R' ^Hn 



R' 



dz 
dr 



R' 



dR 

dr 



R' 



(3.25) 



The origin value of R' (and R') depends on the choice of arbitrary functions, but if r ~ ii there, 

then R' is finite and non-zero, while d'^R/dr^ R" . 

L Jo 

The LT observational relations may be very different from the FLRW ones, especially near the 
maximum in R []93|]. 

Inhomogeneous Models of SNIa Dimming In recent years the LT model has seen quite a 
bit of use in investigations of whether the observed dimming of the supernovae can be explained 
as an effect of cosmic inhomogeneity, rather than invoking a 'dark energy' whose magnitude and 
physical origin are obscure. 

For the case of an observer that is off-centre, [ ]67| ] calculated expressions for the angular vari- 
ation of di, Hq, qq, the source number count, and AT/T, and showed the CMB dipole could be 
explained this way. 



It was first pointed out in Celerier [ ]28| ] that the observed SNIa dimming can be explained by 
inhomogeneity. That paper used a parabolic LT model and showed that a series expansion of ^^(z) 
could easily manifest apparent 'acceleration'. This was generalised to non-parabolic LT models in 



[113]. 



In ^ the authors constructed an LT model that has a low density region (void) at the centre. 



and asymptotically approaches homogeneity. Their functions M & f have the form (2.33) with 



M- 



f2 



H±of3o , 



Aa 
' 2ft, 



1 — tanh 



1 — tanh 



( 2Ar 
/ r-rp 
V 2Ar 



(3.26) 
(3.27) 



so that M goes from ~ H']_Q{aQ — Aa)r^ at the centre to H^^aQr^ at large r, and similarly for /. 
Their 3rd function was fixed via the hyperbolic version of (2.3S), choosing ti to be recombination, 
and setting Ri = a^r where a* is the RW scale factor at recombination. They then calculated 
the redshift and the apparent magnitude for a central observer, and found they could obtain good 
agreement with observations. They also verified that they could retain the observed CMB power 
spectrum. In [Q] and [^] the authors investigated an off-centre observer in two versions of the void 
model. They found a marginal improvement in the fit to the SNIa data is possible. If the observed 
CORE dipole is due to this effect, it requires only a 15 Mpc displacement from centre, but the 
corresponding quadrupole and octopole effects are then too small to match observations. 



In IjnSp it was suggested that inhmogeneous models of supernova dimming have a 'weak 
singularity' at the centre. However, this is merely a conical point in the density profile, and not a 
singularity [^]. Also, [120] showed it is easy to smooth the central density without affecting the 
model much. 

Ref [15 1 considered a selection of LT models, with and without A — a central void model 
with a = 0, a uniform present-day density and a varying Ht, a varying Hr with a = 0, both p 



14 



Modelling Inhomogeneity in the Universe 



and a varying. All the models had small density oscillations imposed to represent smaller scale 
inhomogeneity. It was argued that all of the A = models considered are 'peculiar' . 

In the LT functions were chosen to be 2M = Br^, a = 0, f = r^/ (1 + (cr)^), B was set 



from QmO via ( 2.37| ), and c was adjusted to get the best fit to the SNIa magnitude-redshift data. 



With = 0.2 the fit was better than that of ACDM. 



For observers near the centre of an overdensity, the model of [33] agrees with observations 
for part of the z range. 

[43] also compared two classes of LT models with the SNIa data, calculating x^- One model 
was fixed by Ht = H + Ai^e'^'/'^o and ^^(r) = 2M/{H^qR^) = ^q, where H, AH, tq and Qq 
are constants. The other had a = and 0,^^ varying. They found that varying is very effective at 
fitting the data, but varying Qm is not. The best fit LT model had slightly lower than the ACDM 
model, but including both inhomogeneous expansion and non-zero A did not improve the fit. 



In [ ]22[ ] two models were considered — a local void model with a simultaneous bang time, and 
a 'bubble bubble' model in which the expansion rate Ht is higher locally than far away but the 
present-day density is uniform. Each is a quite specific 2-parameter LT model. They confronted 
their models with SNIa di data, the BAO dilation scale, dy = [d'j^z/Hr]^/^, and the limit on H^ set 
by the age of the oldest stars. From calculations, they concluded that their best-fit hubble-bubble 
model fits the data almost as well as ACDM. 



In 04 5|], void models with 4 or 5 parameters were considered, and it was shown that they can 
provide a good fit observations of the SNIa dimming, the CMB, and the BAO (within 1 a) and a 
comparable to the ACDM model. In [^], it was shown that observations of the kinematic Sunaev- 
Zeldovich effect already limit LT voids to < 1.5 Gpc, and future surveys will either put tighter 



limits on the size or constrain the density and expansion profiles. In []47[], the authors proposed the 
normalised cosmic shear as a test of inhomogeneity. They also found that LT models still provide 
excellent agreement with updated SNIa and BAO data. 

A similar good fit with SNIa observations, i.e. a comparable to that of the ACDM model, 
was found in [^], which considered LT 'bubble' models with decreasing H{r) and constant 
^Imif)- There was no improvement in the fit using a similar model with non-zero A. 

For a summary see [29, 3^, 20] . The important issue here is to highlight the difficulty of sepa- 
rating the effects on the null cone observations of the cosmic equation of state, of source evolution, 
and of cosmic inhomogeneity. Whether or not ~Gpc scale inhomogeneities are discovered, inho- 
mogeneous models have to be taken seriously, firstly because inhomogeneities on many scales do 
exist, and secondly because we should rigorously verify homogeneity (instead of just assuming it), 
and such testing requires using an inhomogeneous model, so that the detection of inhomogeneity 
is a possible outcome. 

Differences between Dimming Modeis Now the arbitrary functions of any given LT model 
determine not only a luminosity or diameter distance relation, ^^(z) or d£i{z), but also a redshift- 
space density relation (t{z). Each chosen model "predicts" a a{z) profile, and this will be important 
in distinguishing models. Though number counts are not very complete or reliable today, the 
situation is likely to improve rapidly with future redshift surveys. In fact, the different types of 
model predict very different a{z). 
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If the LT arbitrary functions are written in the form of a central behaviour plus a variation, 

M = Msr^{l + AM) , f = f2r\l + Af), a = ao + Aa, (3.28) 
then the leading terms may be related to the central cosmological parameters via 

(3.29) 

, ro = M3(sinhr?o-r?o)//2'^'; (3.30) 



/2 = sign(/72(l - 2qo)) , M3 = <7o^o (/2/(^o (1 " 2(?o)))'^' 
/2 > : rjo = arccosh 1 + ■ 



^ ^ /2ffoV/2(l-2go) 



/2 < : 170 = arccos 



«o = to-To 



^ ^ /2gox//2(l-2go) 
90 



ro = M3(r?o-sin??o)/(-/2)^/^ 



(3.31) 
(3.32) 



which defines a 'central RW model'. A model with a pure bang time inhomogeneity, may be 
described by the functions 

M = M3r\ f = f2r\ a = ao + I{e-'-/^ - 1) + Kie-^-^-^^ - I) , (3.33) 
i?o = 0.72, go = 0.22, I = 0.8/Ho, J = 0.5, K = -0.7/Ho, J2 = 0.7 , (3.34) 

and using this, we get good agreement with diiz) from supernova data^. It has become custom- 
ary to compare the measured magnitudes with those expected in the Milne model, i.e. Am = 
m — rn-Miine = 51og(L/LMiine)- The left plot below shows Am{z) against the supernova data of 
Kowalski et al. [|73|], with the blue line for the given LT model, and the red line for the RW 
model with the same central parameters; the middle plot shows the redshift-space density (number 
of sources per steradian per unit redshift interval times mean mass per source), with blue the LT 
model, and red the central RW model; the right plot shows p{to,r) / pcrit,o the density as a multiple 
of the central critical density, on a constant time slice at the present day, against coordinate radius r. 




For a pure mass-geometry-energy inhomogeneity, we can set the LT functions to 

M = Msr\ / = /2r2(l + i?(e-^'/^-l) + G(e^/^^-l)) , a = ao, 

Ho = 0.72 , qo = 0.09 , ^ = 8.2 , F = OA, G= -7.6 , F2 = 0.45 , 

and the following plots show we also get good agreement with the SNIa data^. 



(3.35) 
(3.36) 



For a smoother density profile at the origin, the approach of [ 12C ] may be used. 



Both of these curves have a x that rivals the least squares quadratic fit to the data. 
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There is little difference between the Am{z) curves of these two and many other LT models, 
but the a{z) curves are quite distinct. At present there are many uncertainties connected with 
estimating a; the number count data is not sufficiently complete or accurate, the relation to total 
matter density is not well established, and the evolution of galaxy numbers and masses is an active 
area of research. Still, one may hope this will improve dramatically with the next generation of 
redshift surveys. Other ways of testing models, should be pursued. 

Now, observers are very likely to live on planets, which are highly likely to circle stars in 
galaxies, which have a good chance of being inside clusters within superclusters. In other words, 
many observers are likely to be inside regions of higher density. If the far universe is homogeneous 
(on average), then a model of an observer inside a density peak must have an intervening region of 
lower density to compensate the central overdensity. The varying bang time models that reproduce 
the supernova dimming, also seem to generate a central overdensity quite well. Of course, varia- 
tions in both the bang time and the mass-geometry-energy functions will surely be needed to get 
the best fit to all the data. 

3.3 Determining thie l\1etric of tlie Cosmos 

The Metric of the Cosmos project aims to determine the geometry of our universe from ob- 
servational data. This is an 'inverse problem': given observations such as those described above, 
determine as much as possible about the spacetime metric. In practice, one needs to make assump- 
tions about the cosmic equation of state, etc, but the goal is to reduce them to a minimum. 

An important aim of this project is to determine the degree of homogeneity in the universe. 
The large amounts of cosmological data now flowing in will soon make this a real possibility. In 
order to do this, it is essential to remove the assumption of homogeneity, but since the use of a 
RW model is widespread in cosmological data analysis, many calculations will have to be carefully 
re-worked. Now angular homogeneity — that is isotropy — is easy to check, and does not require 
us to know anything about the PNC. Whatever the variation of observables with z is, it must be the 
same in all directions. But radial homogeneity is not at all easy to verify, since the z-dependence 
of observables depends on several things: the time evolution of the expansion (i.e. the equation of 
state), the source evolution, and whatever radial inhomogeneity is present. So although a general 
treatment requires us to go beyond spherical symmetry, just pinning down the degree of radial 
variation would be a big step forward. 

It was shown in [ ^ ] that any reasonable 'observational' functions doiz) or ^^(z) and a{z) 
can be reproduced by an LT model, and an algorithm for extracting the LT arbitrary functions 
was given. This algorithm was implemented as a numerical procedure and clarified and extended 
significantly in [|7[ |^]. 
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For this project, we now need to invert the equations of the last section; we treat R{z) and 
a{z) as given, and we want to determine /(r), M(r) and a(r). We can use the radial coordinate 
freedom to choose 



dr 



-P{r) , i.e. R' =pW , 



(3.37) 



on the observer's past null cone, so that the solution to (3.9) is 



t = to- 



13 dr 



(3.38) 



In practice, /5 = 1 and t = to — ris the obvious choice, providing there ai^e no shell crossings. Note 
that ( 3.37 ) and ( 3.38D and much of the following only hold for the single null cone with apex (to, 0). 



Putting this in (3.16) gives 



dR 
dr 



■PiW-R), 



and using (3.16), (3.37) and (2.5), we find 



2(3 [dr 



2M _ Aj£ 
R 3 



dR 
dr 



while (3.13) and (3.37) give 



M' = aW 



dz 

dr 



Combining (3.40) with its derivative results in 



(3.39) 



(3.40) 



(3.41) 



Putting ( |3.37D and ( [3.42[ ) in (gJJ) leads to 
dz dr I 1 dR 



dR (3 
■ dr 



/3 dr 



(1 + ^) 



R 



(3 dr 




(3.42) 



(3.43) 



which simplifies, on substituting for M' and R from (3.41) and (3.39), to 



d£ 

dr 



dR 

dr 



(3.44) 



Since the coordinate r is not an observable, we convert all r derivatives to z derivatives using 



dR dR 1 



d'^R d^R 1 dR 1 dip 



dr dz ip ^ dr^ dz^ p"^ dz p^ dz ' 



(3.45) 
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where 



dr 



defines ip. Then ( 3.41 ) and ( 3 .4(1 ) become 

dM 



aW 




2l3(p \dz J 



2(f 



(3.46) 

(3.47) 
(3.48) 



while putting (3.45) in (3.44) and solving for d(f/dz, gives 



d(p 



[l + z) 



+ 



AR 

Az 



Pz_ 

(3 



(3.49) 



where f3z = (3' (p. As stated, the obvious gauge choice is (3 = I, /3z = 0. Equations ( [3.46[ ), ( p.49| ) and 
(3.47) with ( 3.48D constitute the differential equations to be solved for ip{z), r{z), M{z) and W{z). 
Then t{z) and a{z) follow from (p|)-(plll), ( |3l0| ) and ( p8| ). Knowing r(z), M(z), W^(z), and 
a{z) fully determines the LT metric that reproduces the given R{z) and a{z) data. We note that 
( |3.49[ ) is an independent DE, while ( |3.47D and ( |3.46| ) require (^(z). Also ( |3.47| ) with ( |3.48D is a first 
order linear inhomogeneous ODE, for which the formal solution in known. However it is easiest to 
solve all the DEs in parallel as one numerical procedure. 

The initial conditions for these DEs are set at the origin at the time of observation to- The LT 



origin conditions appUcable to these null cone equations were given in [|87|, and are reproduced 
and generalised in the appendix. 

3.4 Apparent Horizon 

In an expanding decelerating model, there is a point on each PNC where the areal radius (i.e. 
dj)) is maximum, dR/dz = 07 We denote this point by R = Rm, z = Zm, and the locus of all such 
points is the apparent horizon (AH) — see [78, 51]. 



But points where dR/dz = make the DEs (3.4S) and (3.47) with (3.48) singular. However, 
in any given LT model is a fixed arbitrary function, so we don't expect any divergence on the 



right of (3.48). Further, it was shown in (3.2C) that 



dR 
d7 







W-R = 



(3.50) 



which open up the possibility that ( [3.48| ) is not actually singular on the AH. Similarly, we don't 
expect dz /dr or d^r / dz^ to be divergent here in a general LT model with co-ordinate choice ( 3.37 ), 
and in fact {^J^ and I ^A^ (with dR/dz = 0) show that 



d^R 



dz^ 



,fR 



R 



(3.51) 



'^This is evident in but first stated explicitly in [ |66| 
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Indeed, (3.50) and (3.51) are exactly what happens at Rm in the FLRW case. So although there 



are no divergencies at Rm, the numerics break down. In 087[ ] this was overcome by doing a series 
solution of the DEs ( 3.46| )-( 3^49| ) in Az = z — Zm, and joining the numerical and series results at 
some z value Zj < Zm — see sections 2.6, 3.3, and appendix B of that paper, and also appendix 
P below. As pointed out in [^], this phenomenon is not merely a cosmological curiosity. At 
this locus, and no other, there is a simple relation between the diameter distance do = R and the 
gravitational mass Mm that is independent of any inhomogeneity between the observer and sources 
at this distance: 



2 Mm — Rr, 



(3.52) 



or R = 2M if A = 0. However, the redshift Zm at which this occurs is not model independent. 
Thus the maximum in R provides a new characterisation of our Cosmos — the cosmic mass. More 
practically, ( 3.52 ) and ( |3.51 ) provide a cross-check on the numerical integration: the M and Lp 
values obtained from the numerical integration must agree with those deduced from the measured 
Rm and am using ( |3.52 ) and ( 3.51| ). This requirement enables systematic errors in the observational 

?oli 



data to be spotted and at least partially corrected, as was done using (3.52) in 
relation ( p. 52 ) generalises to the Lemaitre model, which has non-zero pressure 



In fact, the AH 



Now ( [3.49| ) may alternatively be written as 

dR{l + z) 



6_ 
~dz 



dz ipj3 



'dR{l + z) 




dz (p(3 





R 

1 f dR{l + z) 



-^{l + z)dz , 

R 



(3.53) 
(3.54) 



by ( [3.25[ ), consequently giving 



r{z) = I ipdz 
'o 



dR,^ 

dz 



a 



1-/5 / ^{l + z)dz 
R 



dz . 



(3.55) 



Although this appears to have no singularity at dR/dz = 0, in fact the term in square brackets in 
(3.55) goes to zero, as is evident from ( p. 54) ). 



Some other attempts at solving a version of the cosmological inverse problem [|n|, ^ |1 18p 
got stuck at this locus. See also the discussion of the apparent horizon and the 'critical points' in 



[79]. 



3.5 Numerical Implementation 

Now in reality, the observations do not provide smooth functions R{z) and a-{z), they provide a 
set of discrete measurements of individual sources. In order to proceed, the data must be collected 
into many redshift bins of width 6z, and bin averages calculated. Furthermore, the derivatives, 
dR/dz and d'^R/dz'^ must also be calculated. So, for each of R{z) and a{z) it is necessary to 
fit a smooth function — a polynomial say — to a range of redshift bins, otherwise mild statistical 
variations in R{z) would create wild fluctuations in dR/dz and especially d'^R/dz'^. The degree of 
smoothing is necessarily a compromise between eliminating statistical fluctuations and extracting 
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inhomogeneity. For example, in [90| a quartic polynomial was fitted to 50 bins of width 6z = 0.001 
on either side of the bin in question. 

A second difficulty is that there is no data at the origin itself — where initial conditions for 
the numerical integration are set. The first bin extends from z = to z = 6z, so it provides average 
values at around z = 6z/2. This is resolved by fitting a near-origin series solution of the DEs to the 
first few data bins (see appendix P), and starting the numerical integration further out. 

As explained above, the maximum in R requires a series solution (see appendix In ad- 
dition, it has one undetermined coefficient. Mi. The numerical and maximum-series results are 
joined at some redshift Zj < Zm and this fixes Mi. The numerical integration is resumed at 
2zm — Zj, using the series values there for "initial" numerical values. 

Thus the numerical integration has 4 regions — the origin series, the pre-maximum numerical 
integration, the near maximum series, and the post-maximum numerical integration — which must 
all be properly joined together. 



In 087[ ] the above programme was implemented as a numerical procedure, and tested using 
fake data generated from an LT model. The fake observational data was exact, and contained no 
scatter or errors other than very small numerical errors. Importantly, the numerics successfully 
reproduced the LT arbitrary functions of the various homogeneous and inhomogeneous models 
used to generate the data. This validated the numerical procedure as viable in principle. 



In [90], the effects of statistical and systematic errors in the data were considered. The nu- 



merical program was revised to output uncertainty estimates for each f{z), M{z), and a{z). It 



was shown how to use the data at the maximum of R{z), via (3.52), to detect and correct for sys 



tematic errors in the observational data. Several examples with fake data were given. Lastly, the 



stability of the DEs ( [3.46| )-( |3.49| ) was analysed, and it was shown they are generally stable, except 
for the dM/dz DE which becomes unstable at redshifts larger than Zm- This issue requires further 
attention. 

Application of this method to redshift survey data is under consideration. However, at present, 
available data has a lot of scatter in i{z), 6{z) and n{z), and considerable uncertainty in the source 
properties L{z), D{z) and ^{z) at larger z values. It is particularly likely that studies of the source 
properties at large z have assumed homogeneity, if not a particular FLRW model. A method of 
testing source evolution theories, independently of possible inhomogeneity, was presented in [p^], 
which considered multiple source types and observations at several wavelengths. 
Combining Data witli an Invariant Distance For many purposes, data at the same z are 
grouped together and averaged, and it is assumed the errors cancel out. However, as is well known 
the peculiar velocities of sources create a scatter in z values, especially in clusters, and the ob- 
server's motion creates a dipole, so although redshift can be measured to high accuracy, it is not an 



ideal monotonic measure of distance. According to Walker's and Etherington's argument [ ]119| , |44| ] 
the source area distance = (1 + -z)dD is independent of the observer's motion. Since this dis- 
tance is determined by the geometry of light rays emanating from the emission event, it is also 
independent of the source motion. Therefore, if the data were good enough it would be more 
correct to combine data with the same ds = di/ {\ + z). 

3.6 Clieclcing Homogeneity 

The LT model requires two physical functions of r to be fully specified, so only models that 
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satisfy both conditions (2.31) are homogeneous. As seen above, it is perfectly possible to repro- 
duce one observational function, such as ^^(z), with a variety of inhomogeneous models, so a 
one-function test is not sufficient. Thus a two-function test is imperative for an unambiguous re- 
sult; for example the redshift-space density a{z) can distinguish between models that fit the ^^(z) 
data. Clearly, then, the MoC procedure will provide an important test of homogeneity, when there 
is sufficient observational data of high precision and completeness. If the procedure outputs LT 
functions that are close to the RW form ( 2.31 ) (say within 1 sigma), then this is a strong indica- 
tion of homogeneity. Checking for homogeneity is so important that we should use all available 
tests. Any deep z test of homogeneity will depend on using the correct source evolution functions. 
According to [56], source evolution theories may possibly be tested with detailed multi-colour data. 



3.7 Other Approaches 

Although the above papers are the only ones that are seriously directed at eventually using 
real observational data, [|ll|] did develop a numerical code based on the characteristic initial value 
problem, and [104, ^ also discussed the problem in broad terms. There have been a number 



papers looking at restricted versions of the 'inverse problem' [n9, 34, 120] that only tried to fit the 



diiz), and typically assumed this is identically the ACDM-FLRW curve. Since this only fixed one 



of the LT physical freedoms, the other was fixed by the authors' choice. As already noted, [ 1 1 , 69] 
did not solve the apparent horizon (AH) problem. In [ ]1 18| ] it was mistakenly suggested that it could 
not be solved using inhomogeneous models — see [f?^] for corrections. 



In [g they chose ^^(z) to be that of the ACDM-FLRW model with = 0.3, Vt^ = 0.7, and 
set M = Mqt'^. They considered models with both / = fnwy in which a is not uniform, and a = 0, 
for which / is not uniform. They were able to find a variety of models that solved their inverse 
problem, and some of the varying bang-time models had quite reasonable redshift-space density. 
They did however encounter difficulties at the AH. 

In [0] they assumed diiz) has the ACDM form (i.e. that of an RW model with 0.^ = 0.3, 
JIa = 0.7, fifc = 0), a = constant, and / = Hqt'^ e~'^^^^ , and they succeed in extracting M(r) 
only up to z = 0.4. Their comments below their eq (32) about the inversion method not probing 
the geometry or being unstable, actually originate from the not handling near-parabolic models 
appropriately, and from not identifying the r coordinate freedom. They also seemed unaware of 
earlier work on shell crossings in LT models. 



In [ |120| ] they also assumed the di of a ACDM universe, as well as a = 0, and they used 
the Dyer-Roeder equation to fix the coordinate freedom. They did overcome the AH problem, 
though the details are unclear. Their solution procedure involved integrating outwards from the 
centre and inwards from the AH, so joining the two parts up involved a 'search' through multiple 
integrations to get a matching. They tested different degrees of smoothness at the centre, but 
showed that the results in the outer regions were unaffected. They also investigated the effect of 
a Dyer-Roeder dumpiness parameter that depends on z, and showed that this could reduce the 
amount of inhomogeneity needed to fit observations. 

Evolution of the Redshift Detecting how the redshift of sources evolves with time, may 



become an important method of distinguishing models of SNIa dimming [41 , 82, 1 17, 120]. 
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3.8 General case 



The idea of deducing the metric of the universe from observations was first analysed in the 
classic paper [ ^0| ] by Kristian & Sachs, and followed up in an important review by Ellis et al pO|]. 

There has actauUy been quite a lot of work on this 



Important early ideas appear in [ ]44] , |1 U 
problem JT05| , [T06| , |T07|, [108[ |88|, |89l, |, |, 0, [10|, [l], |5|, g] , especially for the spherically symmetric 
case and its perturbations, though the general case is quite difficult and has not been developed 



very much. See a summary in []57|]. There is much to be done here, especially turning the general 
case into a workable numerical procedure. 



4. The Szekeres Metric 



The Szekeres (S) metric [111, 112] is a very interesting and largely neglected inhomogeneous 
model. Like LT, it is synchronous, comoving, and irrotational, with a dust equation of state. The 
metric is: 



(4.1) 



where e = — 1,0,+1, / = /(r) is an arbitrary function of r, E = E{r,p,q), R = R{t,r) and ' = 
d/dr, and an orthonormal basis for this metric is 



{R' - RE'/E) 



Applying the EFEs, the density and the Ricci and Kretschmann scalars are 

2(M' -ZME' IE) 

up — 



q ■ 



R^{R'- RE'/E) ' 
7^ = 4A + K/9 , 
4 



4Ar . 

+ _[2A + Kp] , 



where 



« - 6^ 



(4.2) 

(4.3) 
(4.4) 
(4.5) 

(4.6) 



The function R{t,r) has exactly the same dynamics and solution (e.g. (2.9)-( 2.16 ) for A = 0) 
as for LT. The function E is given by 



E{r,p,q) = ^ 



S 



S 



(4.7) 



where functions S = S{r), P = P{r) & Q = Q{r) are arbitrary.^ The p-q 2-surfaces and E will 
be interpreted below; in brief the constant time 3-spaces are foliated by 2-surfaces of constant 



The function E is often given in the form 



E{r,p,q)^A{p^ + q^) + 2Bip + 2B2q + C, 



(4.8) 
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coordinate r, which have 2-metrics of spheres, planes or pseudo-spheres, depending on the value 
of e. 

The e = 0, — 1 cases have been ignored until recently [ |62[ [7^ ]. Although quantities like r, M(r) 
etc do not have the same meaning as in spherically symmetric models, curves of constant p & q 
will be called 'radial', 'p-radial' or 'h-radial', and prefixes 'p-' and 'h-' will indicate quasi-planar 
and quasi-pseudospherical quantities. 



See [ [74| ] for a review of its known properties, or [ ]101| ] for an introduction. See also [|60|, |62|] 
for an analysis of the e = +1 and e = 0, — 1 cases. 

Free Functions The S metric has 6 arbitrary functions /, M, a, S, P and Q, which allow 
a rescaling of r plus 5 functions to control the physical inhomogeneity. For the e = case, the 
mapping {S,f,M) — > {S / F, f F'^ , M F'^) for any F{r) does not change the metric, the density or 
the evolution equations. Thus S{r) might as well be set to 1 with F = S. 

Special Cases The S model contains the LT model when e = +1 and S, P, Q are all constant. 
It therefore contains the same special cases, and has geometric possibilites at least as interesting as 



LT. With £" = it reduces to the Ellis metric [38]. It also has a Kantowski-Sachs-type limit, and 



its null limit is a generalisation of the Kinnersley rocket metric [ |55| ] 



No Kiiling Vectors This metric has no Killing vectors [27], but that does not mean it is even 
close to a general inhomogeneous dust solution. It is the r dependence of E that destroys any 
spherical, planar or pseudo-spherical symmetry. Despite the inhomogeneity of the model, and the 



lack of Killing vectors, any surface of constant time t is conformally flat [10]. 

Matching to Vacuum Also, any surface of constant coordinate 'radius' r can be joined to a 

symmetric vacuum metric with spherical, planar or pseudo-spherical symmetry [24, 2^, This 



latter means that the S metric generates a symmetric gravitational field "outside" each and every 
constant r shell. 

Singularities The S model has the same singularities — bang, crunch, shell crossings, shell 



focussings — as discussed for the LT model in §2. 1 . The bang and crunch, where R = 0, still occur 



at t = a and t = a + lirT. Shell crossings however are more complicated, as they occur where 
R' - RE'/E = 0, provided M' - 3M E'/E and e - / are not zero. 

4.1 Riemann Projection 

To understand the metric component {dp'^ + dq^)/E'^, we note that the p-q 2-surfaces can be 
transformed to 2-spheres, planes or pseudo-2-spheres by the Riemann projection: 



-l,E>{): ^ = coth ( - ) cos(0) 



{p 


-P) 




s 


(q 


-Q) 




s 


{p 


-P) 



: coth ( - ) sin((/>) , (4.10) 



-1 , E <0 : - — - — - = tanh ( - ) cos((/)) , 



where A{r), Bi (r), B2{r), and C(r) must obey 

4(AC-B?-B|) = e. (4.9) 

This last is automatically satisfied by ( |4.7[ ), so calculations are easier. Also S, P, Q have a natural interpretation in the 
Riemann projection given next. 
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(q-Q) 



e = : ^ ={l] cos{4>) 



e = +1 : either — = cot ( - ) cos(0) 



or — = tan ( - ) cos((/>) 





s 


ip 


-P) 




s 




-Q) 




s 


{p 


-P) 




s 




-Q) 




s 


{p 


-P) 




s 


(q 


-Q) 



: tanh ( - ) sm(0) , (4.11) 



^)sin(</.), (4.12) 



cot ( - ) sin(</.) , (4.13) 



tan ( - ) sin(0) . (4.14) 



Notice that, with 9 ranging over the whole sphere, each of the spherical transformations (4. 13) 



& (4.14) covers the entire p-q plane. 



In contrast, both of the pseudospherical transformations (4.10) & (4.11), with < < cxo, are 
required to cover the entire p-q plane once; each transformation maps one of the hyperboloid sheets 
to the p-q plane. To distinguish the sheets, we choose 9 to be negative on one and positive on the 
other. Each constant r "shell" seems to be a hyperboloid with two "sheets", but we shall determine 
whether both these sheets are needed, or even allowed. In the planar case, the Riemann projection 
may be viewed as an inversion of the plane in a circle, or as a mapping of a semi-infinite cylinder 
to a plane. These projections are illustrated below. 
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However, the transformations from {p,q) to 



introduce cross terms in the metric, such as 



drdO. Of course, if £" = everywhere, this transformation recovers the LT model. 
4.2 Properties of E 

The function E determines how the coordinates {p, q) map onto the 2-d unit pseudosphere, 
plane or sphere at each value of r. Each 2-surface is multiplied by factor R = R{t, r) that is 
different for each r, and evolves with time. Thus the r-p-q 3-surfaces are constructed out of a 
sequence of 2-dimensional spheres, pseudospheres, or planes that are not arranged symmetrically. 
Obviously, for e = 0, — 1 the area of the constant t & r 2-surfaces could be infinite, but in the e = +1 
case it is ^ttR^. 

In the p-q plane, E has circular symmetry about the point p = P, q = Q, which is different for 
each r. The E = circle 



{p-Pf + {q-Qf 



-eS' 



(4.15) 



has > on the outside, but becomes a point if e = 0, and does not exist if e = +1. The divergence 
of the metric components gpp and gqg as — > has a geometric significance that will be discussed 
below. The E' = locus is also a circle in the p-q plane, which can be written 



p-P P' 

\ 

S S' 



.q-Q Q' 
' S S' 



p'r+m 



r\2 



For e = 0, +1, this locus always exists, and when e = — 1 it only exists if 

iS'f<iP')^ + {Q')\ 



+ e . 



(4.16) 



(4.17) 



with the radius of this circle shrinking to zero as the equality is approached. It can be shown these 
two circles always intersect, if they both exist. 

To see how E' /E affects the metric and the density, we write x = E' /E. Then in the metric 
(P~T|), Qrr is a decreasing function of x provided x > R' /R, while for the density we have 

6M {M'/{3M)-x) 



i?3 {R'/R-x) 



(4.18) 
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so that 

^ dp_ QM {R'/R-M'/{3M)) 

dx~ i?3 {R'/R-xf ' ^^-'^^ 

and if X — > ±00 

Svrp-^. (4.20) 

Therefore at given r and t values, the density varies monotonically with x = E' / E. (The sign of 
the numerator may possibly change as R evolves.) If x can diverge, p approaches a finite, positive 
limit. 

4.3 Spatial 3- & 2-Geometries 

It is apparent from the above that e determines the shape of the constant (t,r) 2-surfaces that 
foliate the spatial sections: 

e = +1 — >• sequence of Riemann spheres 

e = — 1 — > sequence of Riemann hyperboloids (4.21) 

e = — > sequence of Riemann planes . 

This is confirmed by the curvature of the p-q 2-spaces; the orthonormal basis components of the 
Riemann & Ricci tensors and the Kretschmann and Ricci scalars are 

^^(p) (<?)(?) (a) = ' ^^{p){p) = "^^mi) = ;^ ' ' ^^'^^^ 

which also show i? is a scale length for the curvature. In fact, it is quite possible to have the 
three types of foliation in one S model. The original notation [HI, 1 12 | has a continuous function 
instead of e. These 2-surfaces have area 

A = R^jm^, (4.23) 

which is AttR'^ when e = +1, but otherwise is infinite. 

Note that > requires e + / > 0, to keep the metric Lorentzian, and so 

/>0 ^ e = +l,0,-l 

/ = ^ e = +l,0 (4.24) 

-l</<0 ^ e = +l. 

Clearly, the 3-d geometry determined by / may restrict the possible foliations. For example, you 
can't foliate a positively curved space with hyperboloids, but you can foliate a negatively curved 
space with spheres. 

Calculating the orthonormal basis components of the Riemann and Ricci tensors and scalars 
of the r-p-q 3-spaces, we find 

3o _3o _-l {f'/2-fE'/E) 3 _-/ 

^(r)(p)(r)(p) - ^(r)(g)(r)(g) - {R' - RE' / E) ' ^WWWW " ^ ' (^.^b) 
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3p 3p -1 f {f'/2-fE'/E) f\ , 



R V {R'-RE'/E) R 
'2-fE'/Ef f 
i?2 (R' - RE'/EY R^ 
-2 } 2if'/2-fE'/E) f\ 
R V (R'-RE'/E) r) ' 

The flatness condition '^R(a)(b)(c){d) = requires just^ 



, 4_(2_{fr2-JE^ P\ ,427^ 
{R'-RE'/EY ^i?^'' ^ ' 

'n = 4('4/^^^U). (4.28) 



/ = = /' . (4.29) 



This is not possible for e = — 1. When e = 0, ( [4.29| ) would make grr diverge unless R' — RE' /E 



0, which in turn would make p diverge unless M' — 3ME' /E = 0. It will be shown this is only 
possible as an asymptotic limit. 

For the r-p 2-spaces we find 



^R{r){p)(r)(p) - '^R{r){r) - " " 

1 fE,{E'^-E'EjE)-{f'/2-fE'/E)\ ^^^^^ 



R\ [R'-RE'/E) 

dE E -^^ 
dp ^ dq 



where Ep = — , Eq = — . (4.31 ) 



For these surfaces to be flat requires Eg{E'^ - E'Eg/E) - {f'/2 - fE'/E) = 0, and the only solu- 
tion that can hold over an entire 2-surface is /' = = E'. This is because E & E' depend on p, but 
Eg & E'g don't. Note that E" = imphes all of S" = = P' = Q'. Obviously, these surfaces may 
be curved, even when the r-p-q 3-space they foliate is flat. 

RW in Szekeres Coordinates Since the RW metrics can be written in the Szekeres form, it 



is useful to look at the transformations between Szekeres and standard RW coordinates — see [62] 
for a discussion. The k = —1 case allows all three types of foliation, which are compared below in 
a constant t, (p = 0,iT slice. Blue curves are for the e = +1 case, red for e = and green for e = — 1. 
Note that there's distortion, as a negatively curved 2-surface cannot be properly represented on a 
plane — notably orthogonal lines do not look orthogonal. 



^If instead the coordinate dependent condition ^Rabcd = is used, one gets a more complicated result [p2|| 
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4.4 Quasi Spherical Case 



Dipoie The function E describes a dipole distribution [ ]112| , 36, 60] round the 2-sphere at each r 
value, having {E' /E)max = —{E'/E)min located at antipodal points, and £" = on a great circle 
mid way inbetween. From ( p77| ) and ( [4.13 )-(^.14) we find 

S 



E-- 

E': 



1 — COS 9 
S' cos 9 + sin 9{P' cos 4> + Q' sin ( 
1 — cos 9 



(4.32) 
(4.33) 



so the locus E' = 0, 



5'cos(9 + P'sin6'cos(/) + g'sin6'sin(/) = , (4.34) 
is a great circle of the 9-4) sphere. The locations of the extrema of E' /E are found by setting 



d{E'/E)__^ d{E'/E) _^ 



d9 



which give 



tan(/>e = — 



P' 



coscpe = ei 



ei = ±l 



(4.35) 



(4.36) 



29 



Modelling Inhomogeneity in the Universe 



cos 6p 



P' cos <i)e + Q' sin _ ^/W? + W) 



S' 



S' 



£2- 



5' 

£2 



where ei and e2 are independent of e, and the extreme value is 



±1, 



(4.37) 
(4.38) 



E 



-£2 



(4.39) 



extreme 



These two points are symmetrically located relative to the E' = {) circle, forming an equator and 
two poles — a dipole. Naturally, the dipole orientation varies with r. By (4J_) and (16) also 
creates a dipole variation in the ^/g^^ dr metric interval and the density around each constant t-r 
2-sphere. The distance between constant r shells varies with {p,q); RE' /E is the correction to 
the radial separation R' of neighbouring shells, RS' / S is the forward {9 = 0) displacement, and 
RP'/S & RQ'/S are the two sideways displacements {6 = 7r/2, (/> = 0) & (6* = 7r/2, (j) = 7r/2). 



E'/E 



9' 



rr 



P 



max 
min 



mm 
max 



mm 
max 



The interpretation is that the Szekeres 3-spaces are constructed from a sequence of non-concentric 
2-spheres, each having a density distribution that is exactly what's needed to generate a spherical 
field around a new centre. Here we show a section through a set of spheres, the dipole on one 
2-sphere, and a selection of possible r-p surfaces at some moment in time, as well as the dipole on 
a single spherical shell, and some possible r-(p surfaces (of constant t and = 7r/2. 




Regularity The conditions for regular origins, for regular spatial extrema, and for the avoidance 
of shell crossings are similar to those for LT models — see section 22 — except with further 



conditions on the new arbitrary functions S, P and Q. These are laid out in []6Q|], and the case of 
non-zero A is considered in [31]. Near an origin, R—^0, regularity requires 



Mr^R-^ 



f--R^, 5 ~ i?" , P-^R" , Q-^R"" , < n < 1 . 



(4.40) 
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The non-concentric nature of the constant r shells means that shell crossings are more complex 
than in LT models. Adjacent r-shells will first intersect at the point where is minimum, and as 
time goes by the two will intersect on a circle parallel to the E' = great circle — i.e. aligned with 
the dipole. See [50], and appendix ^ for the conditions to avoid them. 

For a regular extremum we require 



f = -l, and M' = f' 



S' = P' = Q' = 



(4.41) 



and the conditions for no shell crossings have to be modified — see [6C]. 

Apparent Horizons According to the standard definition, surface S is trapped if, for any null 
vector field k'^, k^k^^ = 0, we have 



k ;a|2 ^ , 

and the apparent horizon (AH) is the boundary of the trapped region 



;a.\ 



AH 



0. 



(4.42) 



(4.43) 



Now null vectors that are momentarily 'radial', k^ = K{t,r,p,q)(^{R' — RE' /E),jW,0,0) , j 
±1, are also geodesic, k^V^k"- = 0, if 



K' 



K{R' - RE'/E)' 
{R'-RE'/E) 




RE' 



+ 2K\r' 



RE' 



{AAA) 



This together with (2.8) gives 



" R 



R' 



RE' 



R+jW) = ^-^ 



R' 




2M ^ Ai?2 

-Jf + f + ^ +JV1 + / 



(4.45) 



so the expansion of these geodesies is zero when ( |3.52 ) holds and 



—j; that is, for incoming 
rays in an expanding region, or outgoing rays in a collapsing region. Thus [112] found R = 2M is 
the apparent horizon when A = 0. 



The approach in [60] was a bit different, as it was not required that fc'' be geodesic, and it was 
rather determined whether or not null paths were moving to shells of larger areal radius R. It was 
established that, at any given point, the constant r shells are traversed most rapidly by null vectors 
pointing 'radially'. (Null paths that stay radial are not geodesic in general.) It was then found that 
the locus where dR/dr = along a radial null direction (geodesic or not) is not coincident with a 



constant r shell, and is p-q dependent. This was called the 'absolute apparent horizon' in [gOp. Not 
surprisingly, on any constant r shell, the dR/dr = locus is (the history of) a circle aligned with 
the local dipole direction. 

Wormholes We know light can't quite get through the Schwarzschild-Kruskal-Szekeres (SKS) 



'wormhole', and we know that dense LT 'wormholes' are even less traversibile [51]. But, if a dense 



Szekeres wormhole can be bent round round as shown above, so one side is shorter than the other. 



does that make it easier for light to get through on the shorter side? In [60] it was found light still 
can't get through, and ray paths were calculated and plotted for several models. Shown here are 
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some light paths (R) and apparent horizons (A) in a Szekeres 'wormhole', showing fast (f) and 
slow (s) directionf 




-0.25 -0.2 -0.15 -0.1 -0.05 



0.05 0.1 0.15 0.2 0.25 



Now since Szekeres spaces can bend round, this prompts the question of whether we make 
a handle topology by joining up the two sides across a boundary. The Darmois matching condi- 
tions specify how to splice metrics together. (Actually, we don't need the embedding to work, or 
the amount of bending to be sufficient, as long as the Darmois junction conditions are satisfied.) 
However it was found in [ |60| ] that the matching doesn't work, without creating surface layer. This 
result includes the case of spherical vacuum — so the idea of wormhole shortcuts — commonly 
suggested in context of the SKS geometry — is in fact impossible within Szekeres metrics. 

4.5 Quasi-Pseudo-Spherical Case 

When e = — 1, the constant r surfaces are not closed, and the physical and geometric mean- 
ing of R and M have to be re-thought. We lay out some basic properties here, and attempt an 
interpretation later. 

hl-Dipoie Recall that each shell of constant t & r is a two-sheeted hyperboloid of revolution. 
Using ( K.10| ) & ( H-llI ), we can write 

uS 



cosh 6 - 



(4.46) 



E' 



S' cosh d + sinh 6I(P' cos </) + Q' sin ( 



cosh 9 — u 



v = sign(ii^) . 



(4.47) 
(4.48) 



The = circle corresponds to — > ±00, and its neighbourhood represents the asymptotic regions 
of the two sheets of the hyperboloid of (4. IC) and ( [4.1 1[ ). 
The locus ^' = is 



S''cosh6l + P'sinh6'cos(/) + (5'sinh6'sin(/) = , 



(4.49) 
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which only has a solution if (4-. 17) holds. It is a geodesic of the p-q 2-space, and can be pictured as 
the intersection of a plane with a right hyperboloid. 
Writing 



S' cosh 6 + sinh 6{P' cos (j) + Q' sin ( 



the extrema of E' /E sue 



extreme 



^/{S'?-{P'Y-{Q'? 
s 



where €2 = sign (5'). These extrema only exist at finite 9 if 

{S'f>{P'f + {Q'f 



(4.50) 



(4.51) 



(4.52) 



which is the opposite of ( 1-.17| ); so on a given constant r shell, either i?' = exists, or the extrema 
of E' /E exist, but not both. 

It follows from ( |4.50 ) that this extremum is a maximum where E' / E is negative, and a min- 
imum where E' / E is positive. Thus, for each constant r hyperboloid, on the sheet with ES' < 
(i.e. v€2 = —1), E' / E has a positive minimum and goes to +00 as \9\ —>■ 00, while on the sheet 
with ES' > 0, E' /E has a negative maximum and goes to —00. We now specify that < on the 



E <0 sheet. From the foregoing considerations, if ( |4.52| ) holds, then E' /E is the pseudospherical 
equivalent of a dipole, having a negative maximum on one sheet and a positive minimum on the 
other, but diverging in the asymptotic regions of each sheet near E = 0. 



We see in the metric (4J_) that RE' /Eis the correction to the separation R', along the r curves, 
of neighbouring constant r shells, meaning that the hyperboloids are centered differently and are 
"non concentric", as sketched below. In particular RS' / S is the forward displacement, and RP' / S 
& RQ'/S are the two sideways displacements. The shortest h-radial distance is where E' / E is 
maximum. 

■ hyperbola centres 





Regularity Can e = — l regions have R{t,ro) = for some Tq? The derivation of the 'origin' 
conditions ( [4.40 ) does not depend on e, but when e = —1, / ^ is not allowed, since / > 1, so 
'origins' are not possible. 
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The analysis of [62] shows that, to be free of shell crossings, (4.52) must hold. Even then only 
one sheet of the two-sheeted hyperboloid at each r — the one with > (i?V^)max > i^' /E) > 
— oo — can be free of shell crossings. However on that sheet, the conditions are weaker than in LT 
— see appendix As with the LT case, these are obtained by studying the evolution of R' /R. An 
important conclusion is that only one sheet of the Riemann hyperboloid should be used to construct 
regular models, which means not all of the p-q plane is used. 

Regular extrema, R'{t,rm) = 0), are indeed possible. The calculations in |62] lead to 



/ = 1 , and M' = f' 



S' = P' = Q' = 



(4.53) 



but note that the conditions for no shell crossings are more subtle at such a locus. This is not 
an obscure possibility — the k = —1 RW metric in pseudo-spherical coordinates has a spatial 
minimum in R. 

4.6 Quasi-Planar Case 

No dipole For the e = case, we find 



E' 
E' 

'e 



2{S' + 9{P'cos^ + Q'sm(j))) 

S' + 9{P'cos(l) + Q'sm<P) 
S ' 



(4.54) 
(4.55) 
(4.56) 



and though the E = locus has shrunk to the point p = P, q = Q, it still corresponds to the 
asymptotic regions of the plane, 9 = oo. The locus E' = 0, 



S' + P'9cos(l) + Q'9sm(j) = , 
is obviously a geodesic of the p-q 2-space, and it exists provided 

S'y^O and {P' / or Q ' / 0) • 



(4.57) 



(4.58) 



There are no extrema of E' /E, and it its value extends to both ±00. We interpret the above as 
showing that adjacent r-shells are planes tilted relative to each other, with tan(/)o = Q' /P' being 
the direction of maximum tilt, but if E' = Q they are parallel. 

The behaviour found here cannot really be termed a dipole. 
Regularity The question of whether an e = model may have an 'origin', i? — > 0, is a little 



tricky. The origin conditions ( [4.40[ ) require / ^ 0. But if the metric and the 3-curvature is to be 
regular, we expect 



lim Qj. 



lim 

r— >ro 



f 



to be finite and non-zero. Since RE' / (R'E) is not divergent, this implies 



R'r^^r^R 



R' 



b constant. 



(4.59) 



(4.60) 
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while the p-radial distance is 

/Qrr dr ^ r . (4.61 ) 



In other words, R, M and / may only approach zero asymptotically, and the scale of the planar 
foliations becomes ever smaller. This is what happens in the planar fohation of the k = —1 RW 
metric. 

Since the 3-spaces of a completely quasi-planar model consist of planes tilted relative to each 
other, they inevitably intersect somwhere, unless 

S' = P' = Q' = E' = Q , (4.62) 

and the no-shell crossing conditions for / > LT models also hold. This reduces the model to 
planar symmetry — an Ellis model [^]. 

As with the other foliations, regular extrema require no shell crossings and / — > — e, as given 
in appendix But, by (p^-(|42^) and / ^ also requires R' - RE'/E and M' - 



3ME' /E —>■ 0, which is the Kantowski-Sachs type limit for this case. Alternatively, f also 
occurs in the origin requirement i? — > above, which can only be approached asymptotically. 
Quasi-Planar Limit A S model may have both quasi-spherical and quasi-pseudo-spherical 
regions, and the boundary surface between them is a quasi-planar timelike 3-surface. It was verified 



in [ |62| ] that the e = case and projection are suitable limits of both the e = ±1 cases. 



4.7 Phiysicai Discussion of tlie e = -1,0 Cases 
Roie of R 

In the metric ( p~T| ) and in the area integral, A = R^ f l/E'^dpdq, the factor R^ multiplies 
the unit sphere or pseudosphere, and therefore determines the magnitude of the curvature of the 
constant (t,r) surfaces ( [4.22 ). By (125 )-( 4.28), it is also a major factor in the curvature of the 



constant t 3-spaces. Therefore we view it as an "areal factor" or a "curvature scale". However, 
when e < 0, it is not at all like a spherical radius. We note that when e = — 1, there can be no origin, 
but R can have maxima and minima as r varies, while in the e = case, R cannot have extrema, 
and it can only approach zero asymptotically. 



Roie of M In ( [2.5[ ), M looks like a mass in the gravitational potential energy term of the 



evolution equation, while in ( |2.7[ ) M determines the deceleration of R. For e = +1, where the 
surfaces of constant r are spheres enclosing a finite amount of matter, the function M(r) does play 
the role of the gravitational mass contained within a comoving "radius" r. For e < however, R 
is not the spherical radius that is an important part of these ideas in their original form, and M is 
not a total gravitational mass, since the constant t & r surfaces are not closed. Consequently these 
ideas need revising. 

In fact, the impossibility of an "origin" or locus where M and R go to zero when e = — 1 
means that M must have a global minimum, and indeed regular extrema in R and M are possible. 
Therefore the local M value is not independent of its value elsewhere, and integrals of the density 
over a region always have a boundary term, suggesting the value of M (rather than its change 
between two shells) is more than can be associated with any finite part of the mass distribution. 
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In e = models, an asymptotic "origin" is possible, but not required, regular maxima are not 
possible, and regular "minima" are actually asymptotic origins. So, with an asymptotic origin (as 
occurs in the planar foliation of k = —1 RW), the boundary term could be set to zero. 

Nevertheless, the central roles of R and M are confirmed by the fact that the 3 types of Szek- 
eres model can be joined smoothly to vacuum across a constant r surface at which the values of R 
and M must match. The vacuum metric "generated" by the Szekeres dust distribution must have 
spherical, planar, or pseudospherical symmetry, and in each, M is the sole parameter, while R is 
an areal factor. 

We note that, even in the Poisson equation, the gravitational potential does not need to be 
associated with a particular body of matter, and indeed it is not uniquely defined for a given density 
distribution. 

Therefore we find that M is a mass-like factor in the gravitational potential energy. 



Role of / As shown in § |4.3i and as is apparent from the metric (|4.l|), the function / determines 
sign of the curvature of the 3-space t = const, as well as being a factor in its magnitude. In the 
case, e = +1, this 3-space becomes flat (represented in unusual coordinates) when / = 0. In the 
quasi-pseudospherical case, with / = it becomes 'flat' if the signature is made pseudoeuclidean, 
( — h +). In the quasi-planar case, / = is possible as a Kantowski-Sachs type limit. 



As with LT, / appears in the gravitational energy equation (2.5) as twice the total energy per 
unit mass of the matter particles, and we do not need to revise this interpretation. Therefore, this 
variable has the same role as in quasi-spherical and spherically symmetric models. 
Role of E As we have seen, for e = +1, E' /E is the factor that determines the dipole nature of 
the constant r shells, and for e = —1, it is the pseudospherical equivalent of a dipole, except that the 
two sheets of the hyperboloid each contain half the dipole, and only one of them can be free of shell 
crossings. The shell separation (along the r lines) decreases monotonically as E' /E increases. If 
E' = 0, it is uniform, otherwise it is minimum at some location and diverges outwards. For e = 0, 
the effect of E' /E is merely to tilt adjacent shells relative to each other, but only the zero tilt case 
(E' = 0) is free of shell crossings. 

Density Distribution For e = +1 models, the density has a dipole variation around each 
constant r sphere, though the strength an orientation of the dipole varies with r. For e = — 1 models, 
which must have / > 1, it was found that, if /'/ (2/) > M' / (3M) and there are no shell crossings, 
the density is at all times monotonically decreasing with E' /E, but asymptotically approaches a 
finite value as E' / E diverges. Therefore the density distribution on each shell is that of a void, 
but the void centres on successive shells can be at different {p,q) or {6,(j)) positions, in other 
words, the void has a snake-like or wiggly cylinder shape. The minimum density is only zero if 
M' / (3M) = — {E'/E)max- Far from the void, at large 9, the density is asymptotically uniform with 
p & q (i.e. with (/>), but can vary with r. However, where f'/{2f) < M'/{3M) everywhere, an 
initial void in a constant r shell can evolve into an overdensity. 

The no shell crossings conditions imply limits on how far the location of the minimum density 
can be displaced between shells with different r. 

4.8 Applications of the Szekeres Metric 

The Szekeres metric was not used to model cosmological structures until very recently. 
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In 1 12], models of voids next to superclusters were constructed, and it was found that the 
growth of the supercluster is strongly enhanced, relative to the linear perturbation approach. 

In [|l4|], it was found small voids surrounded by large overdensities evolve more slowly than 
isolated voids, while large voids enhance the evolution of adjacent superclusters. 

A swiss-cheese model based on Szekeres inhomogeneities was used in [16] to investigate the 
effect of non-linear inhomogeneities on the CMB. While compensated inhomogeneities have a tiny 
Rees-Sciama effect, the effect of uncompensated inhomogeneities is around ~ 10^'^ and so could 
be responsible for the low multipoles in the CMB. 

The effect of volume averaging was considered in [17] and it was found the results are the 
same as in the LT case, c.f. ||52|]. 



A generalisation of the LT void models for SNIa dimming of §3.2 was given in [[7ip. They 



used a quasi-spherical S model with quite restricted arbitrary functions and few parameters. This 
allowed some angular variation in diiz). It was shown that the model fits the data almost as well 
as the ACDM model, even though the potential of the S model has hardly been explored. 

Shell crossings in certain specific examples of higher dimensional quasi-spherical models 
were considered in [31]. These authors have also investigated generalised quasi-spherical mod- 
els, including collapse and the occurrence of "shell focussing" naked singularities, often in higher 
dimensions and involving non-zero pressure or heat flux. 



5. Conclusion 

The universe is of course very inhomogeneous on many scales. To fully understand how these 
structures evolve, and properly analyse our observations will require the non-linearity of exact 
inhomogeneous metrics. 

Up to now, homogeneity has been assumed, and was key to making progess. In the age of 
precision cosmology, we should thoroughly test this assumption and quantify how good an approx- 
imation it is on each scale. Nearly all data analysis assumes the RW metric. To be sure we avoid 
circular arguments, there is an urgent need to re-do all calculations in a general non-homogeneous 
metric. The methods of inhomogeneous cosmology will be an essential component of this endeav- 
our. 

Lemaitre-Tolman models have produced a wide variety of interesting results, and the investi- 
gations are far from exhausted. 

The Szekeres models have a lot of flexibility, and can be used to model quite complex struc- 
tures — but they have been very little investigated. 

There are plenty of opportunities for good research. 



A. Near Bang and Near Parabolic Series 



Near the bang, where i? = 0, (2.5) is dominated by 2M/R, and the exact solutions (2.9)-(2.16) 
for A = involve the cancellation of nearly identical terms, thus generating large numerical errors. 
Taking our cue from the / = 0, A = solution, (pJOj), i.e. R = {9M{t - 0)^/2)^/^, we write R as 
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a series, 



R = J2Ris\ s = T^/^ = {t-ay/^ , (A.I) 



■t=i 



and put it into ( |2.5| ) in the form 

3RR^ = 6M + 3fR + AR^ . (A.2) 
Solving for each power of s in turn we find 

R = R2S^\ 1 + V \ \ \ \ , (A.3) 

\ 7 63 4 4851 11 7007 2002 J ' ^ ' 

In the case A = 0, this is also the near-parabolic series for small / (and s not small). When A 7^ 0, 
and (2^) is integrated numerically, small / is not problematic. 

It is a good idea to have more terms in the series than the bare minimum, so that there is a 
range where both the series and numerical solutions are accurate, and each calculation provides a 
check on the coding of the other. 



B. Near Origin Series for Observational Relations 

When calculating t, R, z and a for an LT model with given /(r), M(r) and a{r), the origin, 
where all but t and a go to zero, requires special numerical treatment. Therefore it is useful to have 
a series expansion for the null cone quantities in the neighbourhood of the origin. Writing 



R = ^^Rijr'6t^ , 

i=l j=0 



6t = t-to 



we can solve the evolution equation ( |2.5| ) for the coefficients Rij, when j / 0, 

RioA 



Ru 
R21 



l2Ms ARl, 
+ h + 



R' 



R 



12 



\RiQ 2 



20 



M3 

'2Rl, 
RwR2oA\ 1 



+ ■ 



^10 



Ru 



6 

^22 



M4 ^ M3R20 ^ R20X 



2RIq 



Rfo 



(B.I) 



(B.2) 



but the RiQ, or equivalently the origin values of R', R" , R'" , etc, must be found by numerical 
integration, e.g. using ( |Z8| ) & ( p.lS) ). Then using 



i=3 



i=2 



j=0 



i = to + ^ty , 



i=l 



E 

i=l 



R = Y,Rir\ a = Y,Kir\ 



(B.3) 



1=1 



i=2 
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and solving (3.9), (3.11), (3.16) and (3.13) power by power, leads to 

II — — -KlO 5 t2 — — -K20 H 2 ' 

2i?io 2 2 Rf^ 3 J Rii 

Ri = Rio , R2 = R20 — RioRu , 
3M3 3M3 /3M3 



K. 



Rii R^^ \ Rio 

( ^ 

+ -3M3/3 + 4M4/2 + 



+ /5 



2M3M4 , 6M|i?20 



i2 



+ ■ 



10 



2NhRwR2oK + 



ARIqM^A 



(B.4) 
(B.5) 
(B.6) 



1 



(B.7) 



C. Near Origin Series for the Metric of the Cosmos 

Not only do f{z), M{z), a{z) and R{z) all go to zero at the origin, f = 0, z = 0, but more 
importantly we don't actually have any observational data at the origin. Therefore, we fit a series 
solution to the data from the first few data bins. We write the LT arbitrary functions as Taylor series 
in powers of z, 



R = zS — / 

1=1 



r 



^R^z\ a = Y,K^z\ 

1=2 

00 00 

J2^iz\ M = Y,M,z\ f = Y,u 



(C.I) 
(C.2) 



i=l 1=3 i=2 

The coefficients in the R and a series are determined by fitting polynomials to the observational 



data near the origin. The null cone DEs, ( p.46D , ( p.49D , ( p.47[ ), and ( p.48| ), with /? = 1, are then 
solved power by power, from which we find the coefficients of the r, M, and / series. The results 
of a Maple calculation are 



Ri\ 



2 7 



2Ro 



Riz+ { R2 + ] z^ + { R3 + + ] z-" + I Ri + 
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3i?3 K^R 
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24 12 QRi 
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+ — + <^ A"3 + ■ 
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12ft 
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(C.3) 
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Ks IIK2 2K4, 2K2R3 K3R2 2K2RI 29K^ \ 
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+ 



3K2 K3] ARi 

R2 + 2^3 + 2i?4 + } 

^ 10 12 I 3 



K2I AR2 
i?. + 2i?3 + f Y 



Ri J. ] ^'Rj 19K2 K3 
2 9 20i?i 4Ri lORi 3Ri 



[ IIK2 
I 6 



Ro 



{K2 



R3 ^2K2Ra 



ISORj 



K2 



R2 
2i?f 



15 



ml 



3Rf 



(C.5) 



Then a is found from a numerical integration of ( |2.5[ ), using ( |3.10| ), in the form 

a = to — f — T, T= I where = — — ^ — - + {f/z^)-\ — S 



d5 , -2 2(M/z^) , , , 2n 
^, where = ^ ' ' +{f/z^) + —- 
00 <-> -J 



R 

z 



(C.6) 



The accuracy of the series is estimated from the ratio of the last and first terms. If i is the 
maximum acceptable error (comparable with expected numerical error), then the z value where the 
program changes from series to numerical integration is given by 



< L 



Z < 



M3L 

Me 



1/3 



(C.7) 



D. The Near-Maximum Series on the PNC 



Near z = Zm, where the maximum Rm = R{zm) occurs, we can solve the DEs of the PNC by 
writing the LT arbitrary functions as Taylor series in powers of Az = z — z„ 



00 



i=l 



R = Rn, + ^RiAz\ a = Km + ^KiI\z\ r = r^ + ^r,Az\ 

i=2 i=l 

00 00 



(D.I) 
(D.2) 



i=l 



i=l 



The coefficients of the series for f , M, and W are obtained by substituting these series into the 
DEs ( p6| ), ( 0l9| ), ( p7| ), and ^M ), 

again with [3 = 1, and the coefficients in the R and a series 
are found from polynomial fits to the observational data near the maximum in R. Using a Maple 
program, we find 
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<^3 = r4 = < TTfPr- + + 
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M^ = |l-^|:^, (D.8) 
Ml = Ml , i.e. undetermined (D.9) 
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M2 = {i^ + -r^}^-^-^-^, (D.10) 




I KjYi KjYi ( 1 + ) Rm I 3 




(D.11) 
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These are generahsations to A 7^ of the results in [pO[], including a small correction in the expres- 
sion for M4. Note that the coefficients of the r{z) and Lp{z) series are fully determined from the 
data, but in the M{z) and W{z) series, the coefficient Mi remains undetermined. In other words, 
its value is fixed by data elsewhere, not by data at Rm. 
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E. Conditions for No Shell Crossings or Surface Layers 



The following table presents the conditions that will ensure a model has no shell crossings 
or surface layers at any time in its evolution. For e = +1, the first group of conditions are those 
that apply to the LT model, and the second group are the extra conditions needed in the S model. 
The Ellis models [ ^ ] are the e = 0, — 1, equivalent of LT models. Although the no-shell-crossing 
conditions have not been explicitly studied for them, they can be deduced by setting S, P & Q 
constant. 
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M', /', a' 


S', P', Q' 


+1 


>0 


all 


M' > 


5 — 3M 


>o 


f>0 
a' <0 
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S - 2/ 

(no condition where / = 0) 


<0 


T" + a' > 
a' <0 

but not all 3 equalities at once 




= 
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-1 
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<0 
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but not all 3 equalities at once 
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